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Abstract 

We present the results of calculations for Pu and Am performed using an implementation of 
self-consistent relativistic GW method. The key feature of our scheme is to evaluate polariz- 
ability and self-energy in real space and Matsubara's time. We compare our GW results with 
the calculations using local density (LDA) and quasiparticle (QP) approximations and also with 
scalar-relativistic calculations. By comparing our calculated electronic structures with experimen- 
tal data, we highlight the importance of both relativistic effects and effects of self-consistency in 
this GW calculation. 

PACS numbers: 71.27.+a, 71.15.Rf, 71.20.Gj 
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I. INTRODUCTION 



During the last two decades we have been witnessing a surge of activity in many-body- 
theory based methodologies applied to condensed matter physics. Here we are concerned 
with one of them, the Hedin's GW method 1 . This particular approach has not only been 
applied to many different materials but it has also been formally developed with an intent 
to enhance its own applicability or to diagrammatically extend it. 

First applications of GW were of "one-shot" type when one starts with local density ap- 
proximation (LDA) to get one-electron eigenstates and construct the corresponding Green's 
function which is then used as an input to perform only one GW iteration. Commonly, such 
an approach is called Go Wo- It usually improves LDA band gaps in semicoductors^ but has 
an obvious drawback because the absence of self-consistency makes it depending on input 
and not conserving.-^ 

To make the approach independent on the input, the quasi-particle self-consistent GW 
method (QSGW) was introduced a few years ago.-^ In this method the Green's function is 
found self-consistently with approximate Hermitian form of self-energy which is constructed 
to minimize the perturbation while keeping a quasi-particle picture. The approach was 
successfully applied to a wide class of materials including simple metals, semiconductors, 
wide band gap insulators, transition metals, transition metal oxides, magnetic insulators, 
and rare earth compounds. First calculations for actinide metals using this approximation 
and neglecting spin-orbit interaction have also been reported.-^ Recently a scheme based on 
Lowdin's orthogonalization was proposed^ which removes an ambiguity in the construction 
of the effective self-energy in QSGW. The method has also been extended to treat finite 
temperatures^ and to calculate spin wave dispersions^. However, similar to the "one- 
shot" variants of GW, QSGW method is not <3>-derivable^, and, as a consequence, it is not 
conserving. This, for example, results in difficulties to calculate total energy. 

Applications of fully self-consistent GW schemes are not numerous. They have been 
applied for weakly correlated solids^— and for free atoms and molecules^ 1 ^. General con- 
clusion seems to be that for weakly correlated simple solids full self-consistency deteriorates 
spectra as compared to "one-shot" or QSGW approximations but improves total energies. 
For free atoms the conclusion clearly favors fully self-consistent calculations. Based on these 
facts one can expect that in solids the spectra obtained by fully self-consistent GW might 
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be competitive with spectra from QSGW if the corresponding physics is local enough, i.e. 
similar to free atoms. Besides, the fully self-consistent GW is ^-derivable and so it is con- 
serving. Also, it is important to mention the works aimed to enhance the accuracy of GW 
based schemes, their robustness, performance, and convergency issues^ - — . 

Another very active field related to the GW method is its diagrammatic extensions. We 
mention here the approaches which use LDA-based vertex correction^—, the approaches 
which use direct diagrammatic representation for the vertex— ~—, and the approach which 
combines GW and dynamical mean field theory (GW+DMFT)^^. Hedin's equations 
and correspondingly the GW method have also been formally extended to spin-dependent 
interactions^^., to treat the electrons residing in a subspace of the full Hilbert spaced, and 
onto the Keldysh time-loop contour—. 

Very recently the importance of spin-orbit interaction was highlighted for the elements 
with large atomic numbers and it was perturbatively included in "one-shot" GW calculations 
for Hg chalcogenides^. In this work we generalize the GW method to solve equations 
explicitly based on 4-component Dirac's theory, which is important to get meaningful results 
for such elements as actinides. This fact together with uncertainty in respect to what kind 
of self-consistency is better to use for actinides defines the scope of the present work in 
which we apply self-consistent GW method based on Dirac equation to study the electronic 
structure of Plutonium and Americium metals. 

These two metals (especially Pu) have been a subject of intensive studies during last two 
decades. From theoretical point of view the best understanding^ - — was achieved using a 
combination of LDA and dynamical mean field theory^ (DMFT) known as LDA+DMFT 
method. LDA+DMFT calculations have resolved the puzzle of false magnetism in Pu and 
Am metals which appears in density-functional based calculations^— but contradicts with 
the experiments^. 

However, there is a problem with LDA+DMFT type of calculations as the approach is 
not parameter-free and requires the input matrix of on-site Hubbard interactions. On top of 
that there is an uncertainty with double counting correlation effects that are present both in 
LDA and DMFT theories. Therefore there is a significant interest to develop diagramatically 
based approaches such as GW and its extensions that offer the possibility to overcome 
both problems. In respect to Plutonium, our study can be considered as the extension of 
previous work by Chantis et al- who have studied this metal with QSGW without spin- 
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orbit interaction and concluded that correlation effects included in GW make the /-bands 
narrower and decrease the crystal-field splittings as compared to the LDA results. We 
extend the work jg] in three ways: i) include spin-orbit interaction by using Dirac form 
for kinetic energy operator, ii) perform fully self-consistent GW calculation and compare 
it with self-consistent quasi-particle (QP) and local density approximations, and iii) apply 
self-consistent GW method to Am metal. 



II. RELATIVISTIC GW METHOD 



Althought, a truly relativistic treatment of the problem would require the use of rather 
complicated equations of Quantum Electrodynamics, we use a simplified approach. First, 
we neglect relativistic retardation effects in the Coulomb interaction. In this case, Hedin's 
original derivation of his famous system of equations 1 still holds with the only extension 
that all fermionic functions (Green's function and self energy) become 4x4 matrices for 
every pair of space coordinates. Also, in order to perform self-consistent GW calculation we 
need only scalar parts of the bosonic functions (polarizability P and screened interaction W) 
in the space of products of bi-spinors which is similar to the non-relativistic theory with 
collinear spin structures. So, in our method which is described below only the fermionic 
functions (Green's function and self-energy) have bi-spinor arguments. Second, we exclude 
positron states and represent the coordinate dependence of Green's function in terms of 
electron states only 




where k runs over Brillouin zone, is the number of k-points, indexes (A, A') denote the 
electronic Bloch band states, as obtained from relativistic LDA= 9 or Hartree-Fock (HF) 
problem, (a, a') are the bi-spinor arguments, and r is Matsubara's time. Thus, in the 
coordinate-space representation, Green's function is generally 4x4 matrix for every r, r' 
pair. 

Inside the muffin-tin (MT) spheres the Bloch states can conveniently be represented as 
linear combinations of 4-component solutions ip t LE (ar) of radial Dirac equation taken with 
the spherical symmetric part of Hamiltonian inside the sphere 
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^(ar)| t = ^Z t k L V^K), 



(2) 



LE 



where t is the specific atom in the unit cell, L combines all spin-angular quantum numbers, 
and index E differs between ip, ip, and local orbitals. Coefficients Z^£ E ensure the smooth 
mapping between the muffin-tin spheres and the interstitial region as it is standardly done 
in the linear agumented plane wave (LAPW) method. 

In the interstitial region we neglect by relativistic effects, i.e. we assume the small compo- 
nents to be zero and represent the large components of Bloch states as linear combinations 
of two-component spinors 



where G runs over reciprocal lattice vectors; s is spin index, Q is the unit cell volume, 
u s (a) is a two-component spin function, and A^ s are the variational coefficients in the 
LDA eigenvalue problem. We keep the same bi-spinor argument a here with understanding 
that two of four components at every r point in the interstitial region are approximated to 
zero. Such an approximation greatly reduces computational time for GW but is still well 
justified because relativistic effects are mostly confined near the nuclei. We have checked the 
quality of this approximation by performing LDA calculations with and without relativistic 
treatment of the interstitial region, and the differences appear to be very small. 

In our implementation of the GW method we have taken an advantage of the well known 
fact that polarizability and self-energy in Hedin's GW system of equations^ are most eas- 
ily evaluated in (r; r)-representation while the equations for Green's function and screened 
Coulomb interaction are most easily solved in (k; a;/i/)-representation where oj/v denote 
fermionic/bosonic Matsubara's frequencies. So, in our approach we switch from one repre- 
sentation to another using Fast Fourier Transform (FFT) algorithm whenever needed. 

Below we give the most important formulae as they appear in the course of one loop of 
the self-consistency. The expressions (EJ and (j3J) allow us to express G(ar, aV; r) for both 
r and r' being inside the MT spheres as follows (due to the symmetry of the solid we can 
restrict r to be inside the unit cell with R = whereas r' may be inside the unit cell with 
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(^UfcfAr) = E ^(«)G£^^WA(^. ( 4 ) 

EL;E'L' 

Here r is inside of atom t in the central unit cell, r' is inside of atom t' in the unit cell 
R', and the number of different R' is exactly equal to the number of k-points inside the 
Brillouin zone. 

In case when both r and r' are in the interstitial region we have three different repre- 
sentations for Green's function: i) numerical values on regular mesh G(ar, a'r'; r), ii) band 
states representation G k A /(r) which follows from (GQ), and iii) representation in terms of plane 
waves 



G R '(ar,aV;r) = — Ve- ikR ' 

k 

x £ e^ k+G ) r Ms (a)G k G;s , G ,(r)4(a')e-^ k+G ') r '. (5) 

sG;s'G' 

We can easily transform between the representations i) and iii) using FFT while the represen- 
tation ii) is connected to iii) by the formula (|3]). Finally, when one of the arguments (say r) 
is inside the MT space and another one belongs to the interstitial region the representations 
for Green's function are obtained as obvious combinations of the formulae above. 

We begin our GW self-consistent cycle by transforming Green's function from (k, r)- 
representation to the real-space. Then we calculate the polarizability in (r, r)-variables. 
For the r, r' pair of indexes within the MT spheres we have the following expression 



PtLk;t'L'k'( T ) — 

-EE E< M ^W«)i^>)> 

E\L\ E3L3 a 

X E ^tE 1 L 1 ;t'E a L 2 ( T ) E ^tE-iL^t'EiL^P ~ T ) 
E2L2 E4L4 



X 



E^^ 2 (« / )i^(« / )Mi: fc ,), (6) 



where indexes k and k' distinguish bosonic basis functions M\ Lk (product basis functions 
which are scalars in bi-spinor space) with the same angular symmetry and we have omitted 
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argument r of all functions in the integrands. For the MT-interstitial and interstitial- 
interstitial combinations of r, r' we obtain: 

P?L k -Ar) = -££ E<^W«)bW«)> 

E\L\ E 2 L 2 a 

x E G ^x;av(^ t \*' L2;aV (/3 - r). (7) 
= - J2 G ^Ar)G* a ^AP ~ r). (8) 

aa' 

Having calculated the polarizability we transform it to the reciprocal q-space and boson- 
frequency //-representation, which schematically is given as 

FS(r)^Pg{u), (9) 

where indexes i and j refer to the product basis functions. Transformation (Q is performed 
similar to the Green's function which was specified earlier. 

After that, we calculate the screened Coulomb interaction W. It is convenient to divide 
W into the bare Coulomb interaction V and the screening part W: 

wsm = v$ + ws(v). (io) 

In the (q, ^-representation we have to solve the following linear equation system for W: 

k l k I 

Having found it, we switch back from v — representation for W to r-representation and from 
q space to r space. 

Next we find the self-energy. This is subdivided onto three steps: i) we solve an ef- 
fective Hartree-Fock band structure problem which is similar to the familiar Hartree-Fock 
problem but with matrix elements of Hartree and exchange interaction calculated using full 
GW Green's function from the previous iteration. To speed up the process, we calculate 
exchange interaction in real space and then transform it to the reciprocal space and band 
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representation. The solution of the effective Hartree-Fock problem gives us a new exchange 
part of the Green's function G x . In step ii) we calculate the correlated part of the self energy 
S c in (r, r)-representation. Again there are three different cases depending on where r and 
r' belong to: 



^t£iLi;t'£ 2 L 2 ( r ) ~~ ^2 

E3L3 E4L4 kLk'L' 



a 

X G?E 3 L 3 ;t>E 4 LM)W*l,t> k >L>(P ~ t) 

J2^lM'M 2L MWI'' L '), (12) 



X 
a 



E2L2 kL a 

xC 2 ; Q vMC(^-T), (13) 



^ = -G^LAr)WSiP-r), (14) 

In iii) we transform the self-energy back to the band representation in k-space, using the for- 
mulae similar to Green's function. We also transform it from r to Matsubara's w-frequency. 

The last part is to solve Dyson's equation in order to find new correlated part of the 
Green's function G c . We perform this step using band representation in k-space: 

£{<W - Gt(k; u;)E AA „(k; w)}GW(k; u) 

A" 

= G^(k;o;)S c AA ,(k;a;)G^(k;a;). (15) 

This is accompanied by finding new chemical potential //, with total-electron-number- 
conservation condition. Then, we transform Green's function back from (k, u)- to (k, r)- 
representation, and use it to calculate new electronic density and new Hartree potential 
which are needed for the next iteration. This closes our iteration cycle. It is important to 
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mention that we completely avoid convolutions in k-space, which saves a lot of computer 
time as compared to pure k-space implementation. 

We can also perform quasiparticle self-consistent calculations. Different from the QP- 
scGW method by Kotani et al.-), our method is based exclusively on imaginary axis data: 
We approximate frequency dependence of the self-energy by a linear function near zero 
Matsubara's frequency and reduce the problem to the solution of Dyson's equation to one 
matrix diagonalization. Then, as justified in [RefjfJ we neglect by Z-renormalization of 
Green's function and use it as an input for the new self-consistent iteration. 

To get single-particle densities of states (DOS) from the full self-consistent GW approx- 
imation we perform similar linear approximation to the self-energy and compute spectra 
as the final step after the self-consistency is reached. For the low energy behavior of the 
spectral functions this kind of analytical continuation is a lot more stable and produces 
essentially the same DOS as the traditional Pade approximation. 



III. DETAILS OF CALCULATIONS 



Parameters of our calculations are as follows: We use mesh 7 x 7 x 7 in the Brillouin zone. 
Green's function was expanded over Bloch states obtained from LDA based full potential 
LAPW band structures. The number of bands in this expansion varies between 142 and 
168 depending on the k-point in the Brillouin zone. Note that such large number of states 
is only possible only when using real-space based implementation of the GW method while 
using reciprocal space, it is very hard to handle more than 40-50 bands in the LAPW based 
GW method. 

Inside the MT spheres we expand the functions of fermionic type (Green's function and 
self-energy) in spherical harmonics up to l max = 5. Bosonic functions (polarizability and 
interaction) are expanded up to l max = 6. In the interstitial region each function is expanded 
in plane waves. We use more plane waves for bosonic functions (250-300) than for fermionic 
ones. Our full basis size to expand bosonic functions both inside the MT spheres and in the 
interstitials is about 600 depending on the particular k-point. 

All calculations are performed for the temperature 1000K. The LDA calculations use 
exchange-correlation parametrization after Perdew and Wang.— 
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TABLE I: 5/ occupation numbers for 5—Pu (taken at the volume of its <5-phase) obtained within 
scalar relativistic (SR) and fully relativstic (FR) approaches. 

Method 5/5/2 5/ 7/2 5/5/2 + 5/7/2 
LDA, SR 5.17 
LDA, FR 4.15 0.92 5.07 
GW, SR 4.82 
QP, FR 4.26 0.56 4.84 
GW, FR 4.45 0.44 4.89 

TABLE II: 5/ occupation numbers for 5—Pu (taken at the volume of its a-phase) obtained using 
fully relativstic (FR) approach. 

Method 5/5/2 5/7/2 5/5/2 + 5/7/2 
LDA, FR 3.72 1.38 5.10 
GW, FR 4.05 0.81 4.86 



IV. RESULTS 

We first discuss our results obtained by various methods for the number of 5/ electrons, 
n 5 f, as given in Tables HI UTt and IIII1 As follows from experiment^ - — , the 5/ occupation in 
Pu is close to 5, and the corresponding occupation in Am is close to 6. Our scalar-relativistic 
GW result (4.82) is very close to the value 4.85 obtained in the calculation performed by 

TABLE III: 5/ occupation numbers for fcc-Americium obtained using fully relativstic (FR) ap- 
proach. 

Method 5/5/2 5/7/2 5/5/2 + 5/7/2 
LDA, FR 5.36 0.84 6.2 
QP, FR 5.66 0.26 5.92 
GW, FR 5.67 0.27 5.94 
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FIG. 1: (Color online) Total density of states (DOS) of 5-Plutonium as obtained in self consistent 
relativistic calculations. Comparison is made between GW, QP, and LDA approaches. 
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FIG. 2: (Color online) Total density of states (DOS) of Americium as obtained in self consistent 
relativistic calculations. Comparison is made between GW, QP, and LDA approaches. 

Chantis et al.- As it is seen from the calculated data our GW results are consistently less 
than the experimental ones, which may be attributed in part to the fact that we count 5/ 
electrons only inside the MT spheres. In this respect, the LDA results, which are a little 
too large, look less consistent with experiment. There is also a noticeable difference between 
LDA and GW in the separation of n 5 j onto 5/5/2 and 5/7/2 contributions where the GW 
approximation produces more 5/5/2 electrons and less 5/7/2 electrons. An interesting trend 
is seen when one looks at the volume dependence of 5/ counts for Plutonium (Tables [I] and 
ILT1) . Full 5/ occupation is amazingly unchanged but the distribution between 5/ 5 / 2 and 5/7/2 
states changes a lot. 
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We next describe our calculated total densities of states (DOS) (Figures [T]|4|) and partial 
densities of states (PDOS) (Figures [5J1E])- In all plots chemical potential is set to zero. For 
S-Pu (FigJTJ) we notice that the occupied part of the spectrum as obtained using GW, LDA, 
or QP is practically indistinguishable while the unoccupied part is different. Here, the GW 
method spreads the spectrum over a wide energy interval while the LDA produces features 
mostly near the Fermi level. 

The spectrum of Americium (Figj2j) shows no similarity between different methods even 
for the occupied part of the DOS. Here we clearly see the advantage of using the GW method 
which gives the lowest position of the peak at minus 2 eV in much better agreement with the 
experimental valued (minus 2.8 eV) than the LDA or QP approaches do. The unoccupied 
part of the spectrum gets progressively wider when we go from LDA to QP and then to full 
GW calculation. 

Calculated electronic structure of 5-Pu at a reduced volume, corresponding to the volume 
of a-phase (FigJHJ) in general shows broader features than the one obtained for the 8- Pu 
volume. The difference between LDA and GW calculations seems to be reduced. 

The DOS of S-Pu as obtained in scalar-relativistic calculation (FigHJ) differs on a quali- 
tative level from the relativistic result, therefore, it is quite clear that any serious calculation 
for this element should take into account spin-orbit interaction. 

PDOS from all calculations performed in our work tell us that 5/ states in Pu and Am 
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FIG. 3: (Color online) Total density of states (DOS) of Plutonium taken at the volume of its 
o-phase as obtained in self consistent relativistic calculations. Comparison is made between GW 
and LDA approaches. 
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FIG. 4: (Color online) Total density of states (DOS) of (5-Plutonium as obtained in self-consistent 
scalar-relativistic calculations. Comparison is made between GW and LDA approaches. 
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FIG. 5: (Color online) Partial densities of states (PDOS) for Plutonium (taken at the volume of 
its (5-phase) as obtained in self-consistent relativstic GW calculation. 

play a key role in energy region close to the Fermi level. We see the increase in hybridization 
between 5/5/2 and 5/ 7 / 2 states when we go from 5-Pu (Figj5]) to a-Pu (Figj6]), and we see 
practically perfect separation between these states in Americium metal (FigJTJ). 

The difference between QP and self-consistent GW electronic structures becomes more 
clear when we consider the quasiparticle renormalization factor Z (FigJ9l and [T0l) . We cal- 
culate Z factor in band representation according to 



I du ) XX' 
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(16) 
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FIG. 6: (Color online) Partial densities of states (PDOS) for Plutonium (taken at the volume of 
its a-phase) as obtained in self-consistent relativstic GW calculation. 
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FIG. 7: (Color online) Partial densities of states (PDOS) for Americium as obtained in self- 
consistent relativstic GW calculation 



In our case, the indexes (A, A') correspond to the effective Hartree-Fock band structure 
problem where Hartree and exchange interactions are calculated using full GW Green's 
function. The more Z differs from 1 the stronger DOS differs from the effective Hartree- 
Fock band structure which usually has too broad spectral features. On Figj9] and [10] we 
have plotted the diagonal components of Z factor matrices as functions of the band index 
for 80 lowest bands for the k = (0, 0, 0) point of the Brillouin zone. In all cases the position 
of the Fermi level is between band 16 and band 17. Actually there are 6 distinguishable 
bands (5/5/2) below Ef and 8 bands (5/7/2) above Ef which have noticeably smaller Z's 
than the rest of the spectrum. It is also clearly seen that Z's for the /-bands in the QP 
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FIG. 8: (Color online) Partial densities of states (PDOS) for <5-Plutonium as obtained in self- 
consistent scalar-relativistic GW calculation. 
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FIG. 9: (Color online) Band renormalization factor Z as a function of band index for fcc-Plutonium 
(taken at volumes of a- and (5-phases) and for fcc-Americium as obtained in self-consistent rela- 
tivistic GW calculations for k = (0,0,0). 



calculation (0.55-^-0.6) are smaller than those obtained in the self-consistent GW calculation 
(0.65 -T- 0.75). This explains why spectral features in the QP electronic structure are closer 
to the Ef. One can say that in case of Am and Pu the QP approximation looks like being 
overscreened similar to the LDA. 



In conclusion, we have described our implementation of the relativistic self-consistent 
GW method and its application to the electronic structure for Plutonium and Americium 
metals. We have found that the inclusion of relativistic effects in GW is extremely important 
for the proper treatment of the actinides. We also discussed the differrences in spectral 
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FIG. 10: (Color online) Band renormalization factor Z as a function of band index for fcc- 
Plutonium (taken at volume of the <5-phase) and for fcc-Americium as obtained in self-consistent 
relativistic QP calculations for k = (0,0,0). 

functions obtained using the present approach with LDA and quasiparticle self-consistent 
GW approximations. 
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